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OF SMOLUCHOWSKI'S COAGULATION EQUATION 

PETER L. W. MAN*' AND JAMES R. NORMS 5 AND ISMAEL F. BAILLEUL § AND 

MARKUS KRAFT 1 

Abstract. In this paper, two new stochastic algorithms for calculating parametric derivatives 
of the solution to the Smoluchowski coagulation equation are presented. It is assumed that the coag- 
ulation kernel is dependent on these parameters. The new algorithms (called 'Single' and 'Double') 
work by coupling two Marcus-Lushnikov processes in such a way as to reduce the difference between 
their trajectories, thereby significantly reducing the variance of central difference estimators of the 
parametric derivatives. In the numerical results, the algorithms are shown have have a 0(f/N) order 
of convergence as expected, where N is the initial number of particles. It was also found that the 
Single and Double algorithms provide much smaller variances. Furthermore, a method for establish- 
ing 'efficiency' is considered, which takes into account the variances as well as CPU run times, and 
the 'Double' is significantly more 'efficient' compared to the 'Independent' algorithm in most cases. 
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1. Introduction. The simplest of pure coagulation processes puts into play 
chemical species characterised by a single scalar quantity, say their mass, with values 
in a discrete set, say the positive integers. The evolution of the process is modelled 
by a differential equation which gives the time evolution of the concentration fJ-t(x) 
of particles of mass x G N. Given a real- valued function / we shall write (/, /j,*) for 
Ssgn f{ x )^t{ x )- Quantities measured by the experimenter (such as moments) are 
of this form. Smoluchowski 's description of the evolution of /x^ is 



(f,rf) = l Y,{^ x+ y)-^ x )-^y^ K ^y)^ x )^{y)- (1.1) 



d 



At 



The kernel K\{-, •) is a symmetric non- negative function which represents the rate 
at which a pair of particles of masses x and y coagulate to create a particle of mass 
x + y. The term {f(x + y) — f(x) — f(y)} is the change which has occurred in the 
quantity (/, fi^) as a result of this coagulation. The letter A in the kernel stands for a 
e?-dimcnsional parameter. Our aim in this article is to devise a new numerical scheme 
for investigating how the solution ^ to Smoluchowski equation depends on A. We 
shall concentrate on the case of a one dimensional parameter since the same analysis 
applies to the partial derivatives for a multidimensional parameter. 
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There is a large amount of literature concerning the solving of the continuous 
particle sized version of eq. (jl.ip and its many variations such as particle inception, 
surface growth, sintering, and fragmentation [T]-[H]- However, little is devoted to a 
systematic method of sensitivity analysis other than merely simulating the physical 
system in question for various parameter values to measure the change in some quan- 
tity as a result of the parameter change. In this paper, we wish to conduct sensitivity 
analysis by explicitly calculating the parametric derivative of equation (jl.l| . 

There are only a few sources which report on this approach 10]-|13j. One method 
involves a weighted particle method, which assigns to each particle a weight with the 
interpretation of the number of physical particles it represents. This particular method 
[10[ 1 1 2 j considers a finite difference approach where both particle systems (with dif- 
ferent parameters) are simulated together, the only difference being the particles' 
weights. A potentially very powerful method based on the Lagrangian formalism is 
considered in [11]. The idea is to consider an adjoint equation which solves for the 
parametric derivative directly rather than eq. (|1.1[) . This allows the solving of the 
derivative for all values of the parameter simultaneously. 

The aim of this paper is to present two new stochastic algorithms for the calcula- 
tion of parametric derivatives of eq. Ijl.ip with emphasis on variance reduction. These 
algorithms are based on the simple Marcus-Lushnikov process, but we consider how 
two such processes with different parameters can be solved simultaneously in order 
to reduce the estimator variance. These algorithms are presented in section O Their 
mathematical formalisation is detailed in section r2.3.3[ in which we describe how this 
formalism can be used to justify that these algorithms do indeed provide approxima- 
tions of the sensitivity. The quality of these approximations is investigated in section 
where numerical results are analysed. 

2. Central difference estimation of parametric derivatives. As Marcus- 
Lushnikov's process is our main ingredient, let us recall first what it is. Dropping 
the index A, equation Ijl.ip makes it clear that \x t should be seen as a non-negative 
discrete measure on N and (/, fi t ) = J2 x gk /( x )A t *( a; ) as the integral of / against \i t . 
In Marcus-Lushnikov's approach, /i t is approximated by a random finite measure of 
the fornfj] 

1 " 

i=l 



whose dynamics are that of a Markov chain with state space 
Qn ■= In G A4(N) 



^i^^iGN Vt J n=l,2 ) ... I (2.1) 
»=i J 



where TV! (N) is the space of all measures on N. We shall talk of each Xi(t) of S Xi / t ) as 
a particle of the system at time t. Start from [1$ = 5 X4 ; associate to each pair 
(xi,Xj) of distinct particles an exponential random time TSy with parameter ^SiSipll ; 
independent of the other exponential times, and set 

T := min{Ty ; i < j}. 



1 8 x is a Dirac mass at x £ IN. 
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The process fi^ remains constant on the time interval [0, T) and has a jump at time 
T. If T = T pq , set 

= + jr{fixj,+x q - S Xp — S Xq ); 

this operation amounts to removing the particles x p and x q from the system and 
adding the particle The dynamics then starts afresh. 

We shall write /if ' N for the Marcus-Lushnikov process corresponding to the ker- 
nel K\. An obvious way of estimating the sensitivity is to approximate it by the 

(random) ratio (^ +2<i ' N — fi t 2<L ' N )/e. No a priori independence or dependence 

between /if + 2e ' N and /if 2 e ' N is imposed. We are mainly interested in this article 
in producing a stochastic approximation of the sensitivity with a low variance. We 
shall thus try to minimise the variances 

Mi (?) ~ Mt ( x ) \ 



as much as we can for all values of x. For that purpose we shall couple the evolu- 
tion of the two Marcus-Lushnikov processes so as to keep them as close as possible, 

noting that increasing the covariance between /if + 2 e ' N (x) and /if 2 e ' N (x) decreases 

the variance. For notational simplicity, we rename /if ± 2 e ' N (x) as /if' N (x). It is 
perfectly possible to describe both trajectories by the R 2 -valued discrete measure 
(fi^' N , fj,^' N ) G Q 2 N , but to encapsulate the following coupling, we consider the fol- 
lowing approach. 

2.1. Coupling. Here is an example of coupling with framework the unit square 
of the plane. Denote by f(x,y) any probability density on the square, and consider 
the problem of minimizing I := J J Q \x — y\ fix, y) Ax Ay, subject to the condition 
that the two marginals of the probability f(x,y) dxdy on both the x and y axes are 
unifornH. Any measure on the square satisfying this condition is said to realise a 
coupling between the uniform probability on the x-segment [0, 1] and the uniform 
probability on the y-segment [0,1] (regardless of the above optimisation problem). 
The probability dx Ay is such a coupling, but it does not minimise /. This minimum 
is attained for the singular probability on the square with support on the diagonal, 
and uniform on it, signifying maximum correlation between the x and y axes. 

Our framework is more complicated than above as the role of [0, 1] is now played 
by the set of particle approximations, i.e. trajectories ({fh^ }t<£[o.t Qad ]) (f° r some 
final time t cn d) with values in the set Qn of finite measures of the form ^ i S Xi ; but 
the basic idea is the same. The minimisation of / is replaced by the minimisation of 
difference between trajectories as seen in the following paragraphs. 

Denote by and the set of particle^] from fi£ ' N and fJ,^ ,N , respectively. 
(The more particles A t _ and A t + have in common the closer the particle systems are, 
and thus increasing the correlation; the set 

Af:=ATnX+ (2.2) 

2 That is the marginal f(x, y)dy = 1 for each x £ [0, 1], and the marginal f(x',y) dx' = 1 
for each y G [0, 1]. 

3 For a measure = A. y^L-, Sxt , the set of particles is i (5 X1 , . . . , Sx n ) for x± 6 N for all i, thus 
allowing multiple particles with the same size to exist in the set. 




4 



P. L. W. MAN, J. R. NORRIS, I. F. BAILLEUL AND M. KRAFT 



is made up of those particles in commo with corresponding measure nf' N £ Qn 
being the sum of those particles in Xf' . Similarly, we consider 

Xf:=X+\Xf , Xf:=X7\Xf (2.3) 

which are the set of those particles present in Xf but not in XjT (and the other 
way round respectively) — we wish to minimise the numbers of particles in these 
sets. We denote nf' N ,nf' N € Qn as the measures corresponding to Xf,Xf. 
From this definition, we can recover fa + ' N from fa + ' N = nf' N + fi>f' N and similarly 
for /i^' . See that the information held in (Xf,X^~) is the same as that held in 
(X®, Xf, Xf) — thus we seek to describe the full stochastic process by the Revalued 
measure (nf' N , [if' N ,[if' N ) € Q% rather than (fa + ' N , fa' N ) £ Q%- Note that the co- 
agulations for Xf" particles are governed by the rates determined by A ± |e, and that 
certain coagulations, such as between a particle in Xf and a particle in Xf, cannot 
occur in the XT' set because the particle in Xf does not exist in XT' . Furthermore, 
coagulations between a particle in Xf and a particle in Xf cannot occur at all. 

The first version of our algorithm, called Single Coupling Algorithm, tries to 
keep the number of particles from Xf as large as possible, imposing that (as much 
as possible) when two particles, both of which are present in XT' and Xf , are chosen 
to coagulate in one of these systems, they also coagulate in the other. The resulting 
particle must also be present in both XT' and Xf, thus helping to keep Xf large. Of 
course, as the coagulation rates in X t ~~ and Xf differ, we cannot prevent a coagulation 
event of the above kind from happening in only one of the systems; we can however 
minimise the rate at which it happens. 

The Double Coupling Algorithm is a refinement of the previous one in which 
we try to make the creation of particles of Xf and Xf as rare as possible. In addition 
to the above coupling, it considers what happens when a particle from Xf coagulates 
with a particle from Xf^ e (which can only occur in Xf set as mentioned earlier). 
In such an event, the same particle from Xf can be used in a coagulation event with 
a particle from Xf /(S . Out of the three particles from Xf, Xf and Xf, the Xf 
particle contributes size to the other two particles, and is itself removed. More details 
about both couplings are described later. 

One ultimately expects that given enough time, the two systems will behave 
almost independently (i.e. there will be few particles in Xf), but the hope is that 
the divergence in their trajectories is slow enough over the time span of interest. The 
simulation of the sensitivity using two independent Marcus-Lushnikov processes will 
be referred to as the Independent Algorithm; it will be used for comparison with 
the other algorithms. 

Labelling. The usage of the triple measure (fif' N ,fif' N ,[if :N ) 6 Q^ N captures the 
similarities and differences between the fa + ' N and pTT' trajectories. Furthermore, all 
particles are stored one single array, and membership of each particle in one of the 
particle sets Xf , Xf ,Xf is implemented by attaching the particle with a label — 
these being ©, 0, Q respectively. The resulting possible Markov steps are given in 
Tabled 



4 The intersection here is used in the multiset sense - if there are s_ particles of size x G IN in 
X t ~~ and s+ particles also of size x in XT, then the intersection contains min{s_, s+j particles of 
size x. Also, X? in eq. 1 12. 31 1 then has max{s+ — s_ , 0} particles of size x. 
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Table 2.1 
Particle labels and their meaning 



Label 


Meaning 


ffix 

Qx 

Ox 


a real particle (of size x) present only in X t + i. e. present in Xf . 
a real particle (of size x) present only in X t ~ i. e. present in Xf . 
a computational particle (of size x) present in X® i. e. 
a pair of identical real particles, one in X t + and one in X t ~~ . 



Table 2.2 

Possible events described using the labelling notation. 



Type 


Event 


Explanation 


(a) 






[ Qx+y 


if occurs in both Xf and X 4 + 


1 (b) 


Ox + 


Qy ~ 


~* \ Ox © ©y + ©rr+y 


if occurs only in Xf 


(c) 






l_ ©a; © ©y + ©rr+y 


if occurs only in X^ 


2(a) 


®x + 


©»H 


" ©z * ffiz+j/ + ©y+z 


See Double Coupling algorithm ex- 
planation. Only occurs in Double 
Coupling algorithm. 


2(b) 


®x T 


Qy " 


~* ©a:+y © ©y 


This represents a coagulation be- 
tween a pair of particles from 
(X t ®,X t ), so the © particle must 
increase in size and the particle 
becomes a Q (since this particle is 
no longer in the Xf system). 


2(c) 


Qx + 


Qy~ 


~* Qx+y + ©y 


Same logic as 2(b) — this reaction 
can only happen in the X t _ system. 


3(a) 


®x 4- 


©y - 


~* ©x+y 


Particles present in Xf coagulate. 


3(b) 


Qx + 


©y" 


©x-fy 


As in event type 3(a) except for Xf . 



reject 


©x © ©y 


This coagulation cannot occur since 






each of the particles cannot 'see' the 






other. 



Note that there is a certain degeneracy in the state space Q\ — if there exist 
particles j?5 x in both nf' N and nf' N then these two particles can be removed and 
a single particle jt5 x added to nf' N - Note that this cleanup operation is not a 
Markov jump but simply a computational enforcement of the definition of X®, and 
does not affect Xf' at all. 

2.2. Single Coupling system. To be consistent with the above ± notations, 
we shall write K~ for the kernel K x _i e and K + for the kernel K x+ i t . Recall that 
all particles are stored in a single array — we introduce the sets of indices in this 
array which correspond to particles in Xf,Xf,Xf to be I{X?), I[Xf), I{Xf ) C 
{1, . . . , n} respectively, where n is the total number of particleqS We then denote Xj 
to be the size of particle i for i £ {1, . . . , n}. 



5 Thus, I(Xf) U I(Xf) U I{Xf) = {!,..., n}, and their intersections are empty. 
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2.2.1. The idea. The coupling procedure is implemented using a majorant 
kernel. This is a symmetric non-negative function K(-, ■) satisfying K(-, ■) ^ (•, •)■ 
Run both systems X t ~~ and Xf at the same rate, given by K; a coagulation happening 
at that rate is called potential. If a potential coagulation between particles of sizes Xi 
and Xj only happens in X^, perform it with the respective probabilities 

p@ = K +( x » x j) or pe = K~(x i ,x j ) ^ (24) 

KyXijXj) K(xi,Xj) 

otherwise leave the system as it is. This way each system behaves as a Marcus- 
Lushnikov process with the correct rate. 

The coupling itself takes place when the potential coagulation involves a pair of 0- 
particles (where the coagulation can potentially occur in both X t + and X^~ systems). 
In this case, the same uniform random variable on (0, 1) is used to decide whether 
or not we perform the coagulation event in each system. In other cases the potential 
coagulation involves only one system. More explicitly, consider only those pairs (i,j) 
of 0-particles (possibly) involved in the potential coagulation even1@. 

Set 

(i) Kg{xi,Xj) := min{K + (xi,Xj), K~ (xi,Xj)} — rate at which a coagulation 
of the type 0^ + Q Xj -> Qxi+xj occurs, 

(ii) Ag(xi, Xj) := m&x{K + (xi, Xj) — K~ (xi, Xj), 0} — rate at which a coagula- 
tion of the type Q Xi + Q Xj -> Q Xi + Q Xj + ® Xi+Xj occurs, 

(hi) Ag(xi, Xj) := max{K~ (xi, Xj) — K + (xi, Xj), 0} — rate at which a coagula- 
tion of the type Q Xi + & X] — > ® Xi + ® Xj + Q Xi+Xj occurs. 
Figure HO gives a schematic picture of the procedure. 



K° s A+ + A^ 1 



k 


k 




r \ 








O+Q ^ © 


/ 




If K + > K~ 


, + ©^ 


+ 


+ 


Else 


, + ©^ 


+ © 


+ G 



Fig. 2.1. Rate correction for the Single Coupling — generate U ~ (7(0,1) and perform jump 
event according to the given probabilities. 



2.2.2. The Algorithm. Recall the different types of coagulation that can hap- 
pen in the Single Coupling algorithm; they were named 1, 2b, 2c, 3a, and 36. The total 
rate of potential coagulation is defined as 



P ■■= Pi + P2b + P2c + P3a + PSb- (2-6) 



»i.e. , ijei(xf). 
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Algorithm 1: Single Coupling algorithm 



For simplicity of exposition, we suppose particles from Xf,Xf,Xf are all 
stored in one single array, which is indexed by i and j below. 

1 Set t = 0. Set all N initial particles to have labels, 
while t < t en d do 

2 Generate a realisation of the holding time At ~ Exp(p) where p is specified 
in eq. 12.71 and eq. 12.61 and set t <— t + At . 

The following step simultaneously chooses the process k G {1, 2b, 2c, 3a, 3b} and 
the particle pair with indices (i, j) which have the correct labels for the process 
of type k. 

3 Generate an unordered pair of particles with indices (i, j) for potential 
coagulation according to the index distribution 

K(xj, Xj) . . 

2Np ' (2 ' 5) 

where particles i,j do not have opposite signs (i.e. belong to I{Xf) and 
I{Xf) respectively, or the other way round). See section [2",2. 31 for more 
details. 

switch the value of k chosen do 

4 case k = 1 

This case represents the Single Coupling part of the algorithm and the 
following steps JTJT]) exactly identify with Figure |2~T1 

5 Generate random variable U ~ U(0, 1). 
if < K U < K° s then 

6 perform event type la: + Q Xj — * Q Xi +xj- 
else if Kg < KU ^ Kg + Ag + Ag then 

if K+ > K- then 

7 perform event type lb: Q Xt + Q Xj -> Q Xi + Qx 3 + ®x t +x r 
else 

8 perform event type lc: (D Xi + Q X] —> ® Xi + ® X] + Q Xi +x j ■ 
endif 

endif 

9 break. 

For following cases 2b, 2c, 3a and 3b, assume the labels on the particle pair 
match with those as described in steps[T|[T] otherwise swap the indices 

i and j. Let 'w. p.' mean 'with probability'. Recall also the definitions of 

and pq from eq. (|2.4|) . 
10 case k = 2b w. p. p®, perform ® Xi + Q Xj — > (B Xz + Xj + Q Xj - break, 

n case k = 2c w. p. p e , perform Q Xi + Q Xj — *• Q Xi +Xj + ®x r break. 

12 case k = 3a w. p. p e , perform Q Xi + Q Xj — > Q Xi+Xj . break. 

13 case k = 36 w. p. p®, perform ® Xi + ® Xj — > © Xi +x,- break, 
endswitch 



14 If a coagulation occurred, for each particle that has just been involved in 

the coagulation, or newly formed, search for a particle of the same size of the 
'opposite sign'. If there is such a particle (of size x, say), perform a cleanup 
operation: Q x + ® x — > Q x . 

is if there is only one particle left in the system then STOP. 

endwhile 



16 STOP. 
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where the p k represent the majorant rates at which a potential coagulation of type 
k G {1, 26, 2c, 3a, 36} happens: 



Pi == 



1 

2N 



K{xi,x v ) 



Plb ■= 



E/ K( X il X j, 



*°2c := -77 



i£l{X?) 



i,k 
i^I(Xf) 
k£l(X?) 



K(xi,Xk) 



P3a 



1 

2N 



E K( X 3> X 3') 
j,j'&I(Xf) 



P3b 



1 

2N 



E 

k^k' 

k,k'zi{Xf) 



K(x k ,x k > 



(2.7a) 



(2.7b) 



(2.7c) 



Adopting this notation, one can read the details in Algorithm [TJ on page [7] 

2.2.3. Implementation and Complexity. The main implementation issue 
deals with how step Q] of Algorithm [T] is performed. We make the assumption that K 
can be expressed as (for some A): 



K(xi,Xj) = : f a (xj)g a (xj) 



(2.8) 



such a form of K enables an easy performance of step [T] As an example of the 
factorisability condition, the additive kernel K(xi,Xj) :— X(xi +Xj) can be expressed 
as K(xi,xj) — A . Xi + Xj . A, implying that f\{x) — A, g\{x) — x, fiix) = x and 
9i(x) — A. The assumption is not so strict — one need only find a majorant kernel 
with this feature. More details on this majorant kernel factorisation can be found in 
the articles by Eibeck and Wagner [TJ [2] and Kraft and coworkers [3, HJ [14] . 

To see how step [T] is performed, first define C as the set of pairs of distinct indices 
(of particles) such that the pair are not of opposite sign. Thus eq. (|2.5[) can be written 
as: 



K (Xi , Xj 

2Np 



E 



(p,9)ec 



y<™' fa' { x p)9a' i x q) 



E 

a 

E 



E(p, g )6cE tt ' fa'(Xp)g a '(x q ) 
E(p,g)eC fa{Xp)g a {Xq) 



fa{xi)g a (Xjj 

1 

/ fa{Xi) 



(2.9a) 
(2.9b) 



g a {x 3 ) 



E( M )gcEa' /a'feJSa'K) \Ep fa\ x p) E q ; 



(p,«)ec 



g a (x q ) 



(2.9c) 



where E(i,j)ec fa{xi)g a {x ) = J2t fa(xi) J2j ■ (ij)ec 9a(xj) implies the last equality. 
Thus the user must choose a according to the first fraction of eq. (|2.9cp whilst the 
last two fractions are for the generation of the pair of particles. The advantage of 
these methods is two-fold: first we can store the values f a {xi) and g a {xi) in 'binary' 
tree structures which also stores their sums (over i). This allows efficient generation 
from the respective distributions 



Ep fct{ x p) 



and 



J2 Q 9 a (xq) 
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Updating the values in this data structure is efficient. If the number of stochastic 
particles in a binary tree is n, the complexity of updating and generating particle 
operations take O(logn) steps. Furthermore, the generation of a particle pair is 
simple — the factorisation allows one to generate each particle in the pair separately 
meaning that generation of the pair of particles is O(logn) rather than 0(n 2 ). 

In Step [1] of Algorithm [TJ a search of particles for cleanup is required for each 
iteration. This can be achieved by maintaining linked lists of information about where 
particles of certain size and label can be found on the particle ensemble list. In short, 
the Single Coupling algorithm may be faster than the 'Independent' algorithm since 
we need to simulate for one particle ensemble rather than two. On the other hand, 
the cleanup procedure in the Single Coupling requires extra storage of information, 
and computational time to update this information. 

The Single Coupling algorithm is good for the initial prevention of creation of © 
and Q particles, however, as the particle numbers decrease over time, the Single 
Coupling should become less effective. This motivates the Double Coupling procedure 
which is designed to minimise the rate of creation of © and © particles for later times. 

2.3. Double Coupling system. The aim of the Double Coupling algorithm is 
to try to minimise the rate at which particles of type © or © are created; it was briefly 
described in section EOl Figure l^^l presents a pictorial illustration of this coupling. 



Single 



before 



\ 1 



_Q g e 



x y y z x+y y +z 

after 

cancels ify—y ' 

/ \ 

__a © e e 

x y y' z x+y y'+z 



Double (y=y') 



before 



e 


00- 




—>-size 


X 


y 


z x+y y+z 





these cancel 



after 



e e 



x y z x+y y+z 



Fig. 2.2. Pictorial explanation of the Double Coupling algorithm. 

More formally, we 

(i) choose a particle as the common particle for the + and + coag- 
ulations. This is done at the maximum potential rate at which the two reactions can 
occur simultaneously (for a common particle i S /(X®)) 

max < K(xf,Xi), ^ K(x k >,Xi) 

U'£/(Xf) k'el(Xf) 

(ii) choose a © particle k S I(Xf) (for a + coagulation) and a © particle 
j E 7(X®)(for a potential © + + © coagulation) with respective distributions 



K(x k ,Xi) 



and 



K{xj,Xi) 



Sj'ei(xf ) KyZji , Xi) 



(hi) rejection steps are performed to correct the rates according to whether the 
coagulation happens in each of the and X ( ~ systems. 
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Algorithm 2: Double Coupling algorithm - only the part which differs from 
Algorithm [TJ 

The Double Coupling algorithm is almost identical to the Single Coupling 
algorithm, but is modified by replacing cases k = 2b and k = 2c (steps [1] and 
[T]| in Algorithm [T] with a new combined case k = 2 containing the following 
steps (ignore any particle pair already chosen). 

l Choose a particle i £ I{Xf) with the distribution 

T N (+,i) +T N (-,i) = T N (+,i)+T N (-,i) , 

Ei' e 7(xf)[2V(+,i')+TV(-,i')] ' N & 



The first rejection step — the following step is there purely to transform the total 
potential rate of process 2 into J2i>ei(x°) ^v( v i ^ rom Pz- 

2 With probability 

TW(v,») 
2V(+,i)+TV(-,i) 

we continue, else reject by going to Step[5J 

3 Choose a (©, 0) particle pair (j, k) with j 6 I{Xf) and fc £ I(X^) according to 
the respective distributions 

XJ X * Xi) and k S Xk ' Xi) . (2.11) 
T N (+,i) T N (-,i) 



We now have generated a triplet of particles k) e I{Xf) x I(Xf) x I(AT t e ). 
Define the probabilities of the + and + coagulations occurring in X t + , X^ 
respectively as: 

T N (+,i) K + (xj,Xi) T N (-,i) K-(x k ,Xi) 
P(S+q ■= — and Pe+0 := — . (2.12) 

T N (V,i) K(xj > x i ) T N (V,i) K(x k ,Xi) 

The second rejection (steps[2][2|) — this occurs in an almost identical fashion to 
Figure just with different rates. 

4 Generate random variable U ~ £7(0, 1). 
if U< min{p e+Q ,p e+Q } then 

5 perform event type 2a: ® x + Q y + Q z — » ® x + y + Q y +z- 
else if min{p ffi+Q ,p e+0 } < U< max{p ffi+0 ,p e+Q } then 

if P®+G > Pe+Q then 

6 perform event type 2b: (B x + & y — ► ©x+y + Qy 
else 

7 perform event type 2c: Z + Q y — > y + z + © y . 
endif 

endif 



8 Go to Step [T] of the Single Coupling algorithm (Algorithm [T]) 
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2.3.1. The algorithm. The algorithm is the same as the Single Coupling ver- 
sion, except that we merge processes 2b and 2c into a new process 2 whose majo- 
rant rate is f>2 :— p2b + P2c- In describing the algorithm, and given a particular 
particle i G I(A t ), we write Tjv(+,i) for X)j'e/(x e ) K{ x j' ' x i)> an< ^ ^iv( — ,i) f° r 
Sfe'e/(x e ) K( x k' , the maximum of these two quantities is denoted by Tjy(V,i). 
See Algorithm [5] on page [TU] for the description of the Double Coupling algorithm. In 
the next paragraph we directly check that the algorithm produces coagulations with 
the correct rates. 

Double Coupling algorithm rates. In the specification of the state space earlier, 
we recall that pf' N ,pf' N ^p:f ,N are the empirical measures for the Q,©,0 particles 
at time t respectively. Given a particular particle of mass y, the total majorant 
potential rate at which this particle reacts with any © particle is equal to 

f A ,(+,y):=^/?(x,y) M f' JV (x) , (2.13) 



and similarly for Tjv(— ,y) and T/v(V,y), so that these are analogous quantities to 
TN(+,i),TN(—,i),TN(y,i) used in Algorithm [5] Note that T/v(±,y) are functionals 
depending on fj,f^®' N respectively. Also, the majorant rate at which a coagulation 
event of the form ® x + Q y + Z — > ® x+v + Q y + Z occurs i^ 

K°jf{x,y,z) := 

T N (+,y)+T N (-,y) f N (V,y) 

P2 



52[f N ( + ,y)+T N (-,y)}»?' N (y) fo(+,g)+*W(-.y), 
of proccss 2 ^ . lst Rejection step 

Choose particle 

K(x,y) K(z,y) . f K-(z,y)f N (-,y) K+(x,y)f N (+,y) 
■ mm ■ 



T N {+,y) T N (-,y) I K(z,y) T N (V,y) K(x,y) T N (V,y) 



Choose a Choose a Q z Probability of rejecting neither coagulation 

which simplifies to 

K°£ N {x,y,z) = min {r~ (x,y, z),r + (x,y, z)} , (2.14) 

where 

r-(x,y,z) , = K{x,y)K-{y,z) K(y z)K+(x,y) 

T N (+,y) T N (-y) 

Similarly (and dropping the (x, y, z) for convenience), the rate at which only the 
reaction occurs is 

£±]j N ■= [max {r~ , r + } — min {r~ , f + }] l r + <r - 

= max{r~ - r + , 0}, (2.16a) 



7 We use the index D for "Double" ; this distinguishes the quantities to be introduced from the 
similar ones introduced above for the Single Coupling algorithm. 
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and the rate at which only the A t + reaction occurs is 

A£ ,JV :=max{r + -r",0}. (2.16b) 

As verification for the above rate expressions, we note that the rate at which a pair 
of particles (® x , Q y ) coagulates (in X^) is 

E ( k d N ( x 'V' z ) + a d N ( x ,v, z )) v?' N ( z ) = E r+ ( x 'V^ z ) ^f' N ( z ) 

= tf+(a;,y). 

A similar computation is made to check that (0y,Qz) coagulate in Xf at rate 
K-(y,z). 

2.3.2. Implementation and Complexity. Looking at Step [1] of the Single 
Coupling algorithm where the particle pair is chosen, and simultaneously the process 
k, we note that the combined process pi for the Double Coupling is chosen by choos- 
ing either a (0, 0) or a (0, ©) particle pair. Either way, a particle is automatically 
chosen with the correct distribution in equation 1 2 . 1 01 and one of © and is also auto- 
matically chosen with the correct distribution specified in equation 12.111 respectively. 
This only leaves the remaining particle left to be chosen. 

The complexity of this algorithm should be similar to that of the Single Coupling, 
except that the combined process pi requires slightly more work than in the Single 
Coupling. However, the Double Coupling hopefully reduces the number of and 
and therefore would reduce the total rate of reactions. Consequently, there might be 
slightly fewer coagulation events in total. 

2.3.3. Limit coupled processes. Recall the definitions of the measures [if' N , 
pf' N , (if' as given in section 12.11 In the same way as one can prove that the 
Marcus-Lushnikov process converges to the solution of Smoluchowski equation (when 
it is unique jf|, it is reasonable to propose a similar result for the triple of stochastic 
processes (n®' N , [i & ' N , (i e ' N ) . The limiting object (/i®,^ ,/i e ) is a deterministic 
non-negative measure- valued path. Given three bounded functions /, g, h, it satisfies 
the system^ 



^ t (f^?) = l E [f(x + y)-f(x)-f(y)]K° s (x,y)^(x)(if(y) 
J2 fix) [A+(x,y)+As(x,y)] (if [x) (if {y) 



]T f(y)K° D (x,V,z)nf(x)(if(y)(if(z) 

J2 f(v) [A+(x,V,z) + A-(x,y,z)] (if{x)(if{y)(if(z) (2.17a) 



8 See for instance the article [15] of J. Norris or [16] of I. Jeon. 

9 The rates Kp, A^, and T(±,y) are the analogues to K^ N , A^ ,iV and Tjv(±,J/) but with 
dependence on the measures /u? , /uP rather than on nf' N ,nf ' N . See eqs. to 12Hfl l for the 

expressions for K^ N , A^> ,JV and Tjv(±,j/). 
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and 

d_ 

dt 



(9,1$) = ij E 9(x + y)A+(x,y) t if(x)iif(y)+ ^ g(x) Ag(x,y)fjf(x)nf(y) 



2 ^ 

x,y^l 

E 

x,y, z^l 

E 

a:, y, 2:^1 



x,y^X 



[9(x + y)- g{x)] (K° D + A+) (a?, y, z) /if (a) 
A-(a;,y,z) M f(x) /xf(y) /if(*) 



- ^ [ff(a: + »)-fl(x)-5(»)]ir f (a!,tf)Aif(aj)/if(y) 

x,y^l 



(2.17b) 



A similar equation to eq. (|2.17b[) holds for 4f(^h,fif). The reader will get a clear 
insight on the reason why these equations appear by seeing the generator of the 
discrete measure valued Markov chain ^/if' N ,^f' N ,fif' N ) (this will be shown in the 
next subsection 12. 3. 4p . Recall that 



+,N 



,N 



and so it is easily shown that for any bounded functions / and g 



dt 



-y, 
2 ^ 



[f(x + y)- f(x) - f(y)} K+(x, y) rf{x) fifiv) 



and 



^(ff, Mt ) = \ E ^ x + y)~ 9{x) - g(y)] K (x, y) n t (x) /i t (y) 



This is in accordance with the fact that are Marcus-Lushnikov processes with 

rates K + /~, therefore their limits are solutions to Smoluchowski equation with the 
corresponding rate (under certain conditions). This implies that their difference con- 
verges to the difference of the two solutions, this being true independent of the cou- 
pling. 

2.3.4. Generator. This section gives a description of the generator of the Markov 
chain corresponding to the Double Coupling Algorithm. The stochastic jumps are 
described by the following elementary operations on measures corresponding to the 
jumps indicated in Figure l2~2l We adopt the notations /t for a generic element of 
and x,y,y',z for integer masses. 



J(i(^,y,y' 
Ji c (^,y,y' 



J^,{^,x,y,z 

J3a(V,%,y 

Jil(v,x,y 



f* 1 N 

» + M 






, 6y+y> 


- Sy - Sy, 


, 




Sy + y' 


, Sy - 


Sy, 


. 5 V + 


Sy, 


Sy + Sy, 


, Sy - 


Sy, 






~0 X ~\~ X -\-y 


, Sy 




, S z 


f Sy + z 


—S x ~\~ $x-\-y 


, Sy 




> s v 




Sy 


, Sy 




, -s z 


+ Sy + Z 





, 




, $x+y 


-S x — S 


$x-\-y 3 X Sy 


, 




, 
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The introduction of the following notation clarifies the description of the generator of 
the Markov chain corresponding to the double coupling algorithm. For any 7 6 Qn-, 
define the rescaled counting measure 7 6 Qn on ordered pairs of masses of distinct 
particles as 

7(A x A') := ~/(A)f(A') — H A') , A,4'cN. 

Set also for any = (0i, 02, <fe) and /x = (/Zi,/z 2 , /U3) £ 

(0,M> : = ((01, Ml), (02, M2>, (03,M3» • 

For any measure /1 := (/x®, /i®) € Q"^, with the corresponding seis of particles 
being (A" e , A" , X e ), we have 

((<h,<h,<M,G t N 0i)) := ((<h,G?' @ W) > (^2, af' (M)) , (03, Sf' e (M))) 
where 

(0i,af ,9 V)) 

= ^ E <M?y + y')A^(y,y)M (2/,j/)+ I] 0i(y) AJ (y, y V 3 2/0 

+ E [<M* + v) (n + A+jfx.^l/W^f^ 9 ^) 

+ ^ 0i(y) A^(a:, y, z) n®(x) fi G (y) ^ e (z) 

+ \ E [0i(^ + ^')-0iW-0i(^)] ^ + (^^')^ e (^,^) (2.18a) 

and 

(0 2 ,af' (M)) 

= ^ E [02(2/ + y , )-0 2 (y)-0 2 (y')]^5(2/ I y')^ (2/,y') 
- E 2 (2/) [A+(i/,i/ / ) + As(i/,tf')] Jl e (y,y') 
E 02(jy)A'°(x,y,z)^(x) M Q ( 2/ ) Ai e (z) 

E 02(y) [A+( 2 ;,y,z) + A D (x, 1 y,z)] »® (x) »° (y) » e (z) , (2.18b) 
with , C/ t ,e (/i)> being defined analogously to eq. (|2.18a|) . 
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3. Numerical Results. The results presented consider two kernels: the additive 
kernel K{x,y) — \(x + y) and a kernel that is used in modelling soot formation in a 
free molecular regime (thus we shall call it the 'Soot Kernel' 



K{x,y) 









(=■ 











The reference value of A for the additive kernel will be 1 and for the soot kernel 2.1. 
We shall always take as initial condition for the Marcus-Lushnikov process N particles 
with mass equal to 1. Throughout this section we shall denote by N the initial number 
of particles in each system (which is the same) , by A the above reference value of the 
parameter, whose perturbation will be denoted by e (i. e. the systems are governed 
by the parameter values A ± he), and by L the number of simulations with the same 
initial conditions. The remaining notation is given below. 

(i) t — time of evolution of the particle system 

(ii) t run = time taken to run the algorithms (CPU time). 



(hi) The estimate of (/, /z*) given by the I 
where / is a suitable test function. 



th 



simulation is denoted by F, 



(iv) The estimate of -§^{f, Mt) given by L simulations is denoted by F ; it is 
equal to xEf=i F / A) - 



3.1. Some initial plots. Figure [37T] shows what the derivative of the parametric 
solution of Ht(x) looks like for the Soot kernel for two different evolution times t. 
Figure [22] shows similar quantities, but for the Additive kernel; it is in good agreement 
with the analytic solution given in [2]. 



2 "> 
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(a) t = 0.5 



(b) t = 4.0 



Fig. 3.1. Derivative -^^(x) versus particle size x, for Soot kernel using the Double Coupling 
algorithm, A = 2.1, e = 0.03, N = 10 7 , L = 300. Confidence intervals have been omitted since they 
negligible. 



3.2. Convergence study. There are two sources of systematic error in using 
the central difference estimator — one due to using a non-zero value of e and the 
other due to assuming a finite particle system. It is the latter we investigate — here 
we estimate the order of convergence of the systematic error as iV varies. The value 

10 This kernel is studied extensively in [3] and used in 1171 1181 [T9l . 
1:L Note that taking f(y) = l y=x gives -^p,^(x). 
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3, *~ 

S 7 " 



12 3 4 

log(size) 




12 3 4 

log(size) 



(a) t = 0.5 



(b) t = 2.0 



Fig. 3.2. Derivative gj/llj (a;) versus particle size x, for additive kernel using the Double Cou- 
pling algorithm, A = 1, e = 0.06, N = 10 6 , L = 300. The line joins the points of the derivative of 
the analytic solution of eq. HI. Il l with monodisperse initial conditions, as given in Confidence 
intervals have been omitted since they are visually negligible. 



of e will be fixed here. 



We define the systematic error due to N as the difference between the expected 
central difference (for finite particle number TV) and the analytic central difference. 

e sys (iV; A, e, t) = EF (A ' e) (TV; t) - f^\t) (3.1) 

However, as we nearly always do not know the analytic central difference f^' 6 ', we 
estimate it using 

F (X ' e) (iV larg0 ;t) (3.2) 

for very large Marge- Also, EF (A,e) (7V;t) i s estimated by F {X ' e) (N; t). Now we set 

test function f(y) = li y =i\ again, to ensure that F *' e \N;i) is an estimate of the 

number density for particle size i. Therefore we rename F " '(JV;t) as F (N;t) 
for number density estimate at particle size i € N; we adopt analogous notation 
y>(A,e,») anc j 6sys (7V; A, e, t, i) for this particular choice of test function. Our metric 
for considering convergence in TV is simply the absolute estimated systematic error, 
summed over chosen evolution times (t/s^—jO and summed over particle sizes i: 

T 

ctot = ^2^2 |e sys (iV; X,e,t k ,i)\ (3.3) 

k=l ieN 

Figure l3~3l shows what we expect — the ctot ~ j?- We obtain similar plots for different 
values of e. 

3.3. Statistical error. This quantity is defined as 

e stat (TV; t, A, e) = F (X ' e) (iV; t) - EF (A ' e) (N; t). (3.4) 

12 We take the t k to be (0.5, 1.0, . . . , 7.0) 
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(a) Additive kernel, X = 1, e = 0.06 (b) Soot kernel, A = 2.1, e = 0.03 



Fig. 3.3. log ctot versus log AT, Af = 25 x 2' /or i = 0, . 
Confidence intervals given are for the Independent case only. 



, 7, N X L = 10 8 and < i < 3. 



This is a signed measure which associates to each i 1 the number e s t a t(-^j ^ ^> e ; 
We consider in this paragraph how it behaves according to the different algorithms 
and kernels. The variance of each estimator F\ can be estimated for each i 1 by 



,W ._ 



^X>(z)-F(z)) 2 . 

1=1 



(3.5) 



This implies that the asymptotic 100(1 — a)% confidence interval for MF^'^ (N; t; i) 
is: 



F (X - e \N;t;i)~z a/2 ' '' 




=;(A,e) 



, F v ' ; (A/;t;i) + z Q/2 1 



L 



(3.6) 



where z Q/ /2 is the upper a/2 point of the standard normal distribution. Hence, 



/ w 

P(|e s tat(iV;i,A,e)| < z a/2 \l wl 



(3.7) 



For this paper, we set a = 0.05 i. e. we consider 0.95% confidence intervals. Also, 
consider the sum e to taistat of the single-simulation variances over the particle sizes: 



^totalstat : _ 



J2< 



(3.8) 



We wish to see how this quantity behaves with N. Figure [3741 demonstrates that for 
all three algorithms, the total variance etotaistat behaves as since the slopes of the 
fitted lines are approximately —1. More importantly, the intercept for the Double 
algorithm is lower than the Single case, and much lower than that of the Independent 
case indicating that etotaistat is much smaller for the Double and Single cases than for 
the Independent — etotaistat for the Independent case at one point is approximately 20 
times larger than that for the Double case at t — 1.0. One interesting observation is 
that the difference in the intercepts (of the fitted lines) between the Double and Indep 
decreases, showing that the benefits of smaller statistical error in the Double case 
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become less pronounced as t increases. The same happens for the Single algorithm, 
but at a faster rate, showing that the X t ~ and systems diverge from each other, 
but faster for the Single algorithm than for the Double algorithm. 
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Fig. 3.4. log e to taistat versus log TV over different values oft for all the algorithms for the 
Additive kernel, where N = 2 10 , 2 12 , 2 14 , 2 16 , 2 18 , N X L = 2 26 and e = 0.06. 



3.4. Efficiency. In this subsection, we will discuss which algorithm is 'best'. A 
quantification of this quality needs to be defined which will take into account the 
accuracy and the run time of each algorithm. One such measure can be described as 
the run time needed to achieve a certain fixed statistical error (assuming N and e are 
fixed). Before defining it, we shall 

(i) Let t Tun (t) be the CPU time actually taken to perform L run simulations. 

(ii) Set efixod = 6tot ^ l3tat ; this "total standard error" will be artificially fixed. 

(iii) Let L cst (t) be the estimated number of simulations required to acquire the 
fixed efi xo d(i) value and t est (t) be the CPU time required to perform L cst (t) runs. 



1 F 



implies that L est (t) = 



we then have 



The condition e fixed (t) = x; sU 

that t cst (t) = ['°'^ £ cst (t). We shall use the following quantity to compare different 
algorithms. 



Inefficiency algorithm := 



t 



algorithm 



est 

y-Doublc 



est 



(3.9) 
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we call it the Inefficiency with respect to the Double Coupling algorithm. 

Hence, if an algorithm has an Inefficiency of more than unity, the algorithm does not 
perform as well as the Double Coupling algorithm. Figure RTijl plots these inefficiencies. 
One can see that the Independent algorithm has large inefficiencies for small t — 
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Fig. 3.5. Additive kernel — Inefficiency relative to Double for Indep and Single algorithms, as 
a function oft for different values of e and N. The Independent algorithm is represented by circles, 
Single by triangles and the Double threshold by the horizontal line. 



this is due to the vastly smaller statistical errors of the Double and Single Coupling 
algorithms, as well as all three central difference algorithms taking comparable times 
to run. In fact, the Double and Single algorithms are generally quicker for smaller 
e since the Independent algorithm requires two simulations to generate a derivative 
estimate, as well as the fact that the Single and Double algorithms have fewer and 
© particles to deal with. The inefficiencies of the Single algorithm lie between 1.0 and 
2.0 (note that Figure [3T5l uses log scales) indicating that the Double algorithm has a 
significant improvement over the Single in terms of accuracy. Also, the inefficiencies 
decrease with e since larger e implies that © and particles are created, thus meaning 
that both coupling algorithms have larger CPU run times (the run times increase 
almost linearly with e for the Single and Double, whereas they are almost constant 
with respect to e for the Independent algorithm). Also, the ratio of variances for 
the Independent and Single algorithms relative to the Double algorithm decreases for 
larger e, due to the increasing similarity between the Independent algorithm and the 
coupling algorithms, thus decreasing the inefficiencies. Note however that the case 
e = 0.2 is unlikely to be computationally useful as it amounts to a 20% change in the 
parameter, which enters multiplicatively into the kernel. 
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Furthermore, one notices that the larger value of N — 10 5 results in larger in- 
efficiencies for the Independent and Single algorithms for the Additive kernel. This 
appears to be because as N increases, the capacity for cancellations is larger, implying 
that the Double is more accurate (and faster) than one expects. Also, as t increases, 
we find that the number of particles for all algorithms decreases dramatically, and so 
there appears to be little difference in accuracy between all three algorithms for larger 
t. This is to be expected since and [l[ will become increasing dissimilar, thus in- 
dicating smaller covariances, and so the variances for the Double and Single Coupling 
algorithms will be similar to those for the Independent algorithm, this loss of efficiency 
of the coupling algorithms is ultimately unavoidable in this class of algorithms, and 
one can only hope to minimise this decrease. 

It is important to realise that this analysis does not take into account the system- 
atic error due to e or N since the Inefficiency metric only uses estimated variances. 
A related problem with the analysis is that the number of particles for t S [3.5, 7.0] 
becomes quite smal£f|, and therefore the systematic errors and estimated variances 
are not very reliable. 

4. Conclusions. In this paper, two new stochastic algorithms were described 
which solve for parametric derivatives of the solution to the discrete Smoluchowski's 
coagulation equation. These algorithms consider two Marcus-Lushnikov processes 
which are coupled together in order to reduce the difference in their trajectories. The 
hope was that this would significantly reduce the variance of the central difference 
estimators of the parametric derivatives. In the numerical results section, we first 
validated the fact that the order of convergence for these algorithms is indeed 0(1/N). 
Furthermore, it was shown from the statistical error plots that the accuracy is order 
of magnitudes better than that of the worst case (the Independent algorithm) , at least 
for larger N and smaller e. Subsequently, we considered a method of comparing the 
algorithms which considers both the variances of the derivative estimators as well as 
the CPU run times. It was shown that the Double algorithm is mostly more 'efficient' 
than Single over variations in e and t, whilst being significantly more 'efficient' than 
the Independent algorithm for small t, large N and small e, though some of this 
advantage is lost for larger t and e. 
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